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Simple practical expressions are put forward, which allow to estimate thermodynamic properties 
of Yukawa fluids in a wide range of coupling, up to the fluid-solid phase transition. These expres¬ 
sions demonstrate excellent agreement with the available results from numerical simulations. The 
approach provides simple and accurate tool to estimate thermodynamic properties of Yukawa fluids 
and related systems in a broad range of parameters. 
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I. INTRODUCTION 

Studies of static and dynamical properties of Yukawa 
systems constitute an import interdisciplinary topic with 
applications to strongly coupled plasmas, dusty (com¬ 
plex) plasmas, and colloidal dispersions. An idealized 
model deals with point-like charged particles immersed 
in a neutralizing medium and interacting via the pair¬ 
wise potential of Yukawa (Debye-Hiickel) type, 

V{r) = (QV»’)exp(-r/AD), (1) 

where Q is the particle charge, Ad is the Debye screening 
length associated with the neutralizing medium, and r is 
the distance between a pair of particles. Clearly, such an 
idealization oversimplifies considerably the actual rather 
complex interactions between the particles in real sys¬ 
tems, in particular in dusty plasmas and colloidal suspen¬ 
sions [iHdj. Nevertheless, many experimentally observed 
trends can be reproduced by this simple consideration, 
at least qualitatively. Hence, it can be considered as a 
basis for constructing more realistic models. 

Thermodynamic properties of Yukawa systems are rel¬ 
atively well investigated using various computational and 
analytical techniques. Some relevant examples include 
Monte Carlo (MC) and molecular dynamics (MD) nu¬ 
merical simulation s. ITI-Il^ as well as integral equation 
theoretical studies |lll - [l4]| . Semi-empirical fitting formu¬ 
las [iM3 and simplistic approaches to estimate 

thermodynamics of Yukawa systems have been also dis¬ 
cussed in the literature. Their accuracy is in most cases 
not better than moderate. 

The purpose of the present paper is to put forward 
practical approach to evaluate thermodynamic properties 
of Yukawa fluids across coupling regimes. This approach 
is based on simple phenomenological arguments, which 
are likely applicable to a wide class of soft repulsive in¬ 
teractions. In particular, expressions for the internal en¬ 
ergy, pressure, and isothermal compressibility modulus 
are derived. When compared with the “exact” results 
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from numerical simulations, the approach demonstrates 
an impressive accuracy. Hence, it represents very con¬ 
venient tool to estimate thermodynamics of Yukawa and 
other related fluids. Among expected applications, wave- 
related phenomena in strongly coupled complex (dusty) 
plasmas can be particularly mentioned. 


II. MODEL 

We consider N particles contained in the (three dimen¬ 
sional) volume V and interacting via the Yukawa poten¬ 
tial O- The number of particles and the volume are very 
large, while the number density n-p = N/V is finite, so 
that finite-size and surface effects are not important. The 
system of repelling particles is stabilized by the presence 
of the neutralizing medium. The effect of this neutral¬ 
izing medium on the thermodynamic properties of the 
system can be trivially evaluated, as discussed in the Ap¬ 
pendix |B] We focus therefore on the contribution coming 
from particle-particle interactions and resulting correla¬ 
tions. Thus we effectively consider a single component 
Yukawa system. 

The system is characterized by two dimensionless pa¬ 
rameters. The first is the coupling parameter, T = 
Q^/oT, where a = (S/dyrrip)^/^ is the Wigner-Seitz ra¬ 
dius and T is the temperature (in energy units). This 
is roughly the ratio of the (Coulomb) interaction energy 
between neighboring particles to the kinetic energy. The 
second is the screening parameter k = a/ Ad, which is 
roughly the ratio of the interparticle separation to the 
screening length. 

The main thermodynamic quantities considered here 
are the internal energy C/, Helmholtz free energy F, 
pressure P, and the isothermal compressibility modu¬ 
lus pL = T~^{dP/dnp)T- In reduced units these are 
u = U/NT, f = F/NT, and p = PV/NT (the ratio 
Z = PV/NT is also known as the compressibility factor), 
respectively. Except explicitly specified, we consider only 
the contribution from the particle-particle correlations. 
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III. DERIVATION OF ENERGY AND 
PRESSURE 

The reduced excess (over that of non-interacting par¬ 
ticles) energy can be divided into the static and thermal 
components 


— '^st “t” ^th- 


( 2 ) 


The static contribution corresponds to the value of inter¬ 
nal energy when the particles are frozen in some regular 
structure and the thermal corrections arise due the devi¬ 
ations of the particles from these fixed position (conse¬ 
quence of the thermal motion). Of course, such a division 
is only meaningful when the regular structure is specified. 
For crystals, the corresponding lattice sum is a relevant 
choice fort Ust- For fluids, it is convenient to link Ust with 
the energy obtained with the Percus-Yevick (PY) radial 
distribution function of hard spheres in the unphysical 
limit 77 = 1, where rj is the hard sphere packing frac¬ 
tion [^, . For Yukawa system this is equivalent to the 

result of the ion sphere model (ISM), where each particle 
is placed in the center of the charge neutral Wigner-Seitz 
spherical cell and the energy is then calculated from sim¬ 
ple electrostatic consideration M- With this choice, the 
static component of the internal energy of the single com¬ 
ponent Yukawa fluids becomes 


Ust = Mf{k)r 


k{k + 1)F 

(k + 1) -I- (k — l)e^'' ’ 


(3) 


where Mf has been termed the fluid Madelung con¬ 
stant [2l|. 

As Rosenfeld and Tarazona [13, [HI first pointed out, 
the thermal component of the internal energy exhibits 
quasi-universal behavior for a wide class of soft repulsive 
potentials, including the Yukawa case. This is illustrated 
in Fig. [1] where the dependence of Uth on F/Fm is plot¬ 
ted for a number of screening parameters k < 5. Here 
Fm is the value of the coupling parameter at the solid- 
fluid phase transition (melting). To produce this plot, 
the numerical data on Uex(K,r) and rm(K) tabulated in 
Refs. [3, [13 have been used (the contribution of the neu¬ 
tralizing medium has been subtracted). It is seen that 
numerical data have a tendency to group around a single 
quasi-universal curve. Reasonably accurate fits (shown 
by the curves) can be obtained by using the following 
functional form 


uth = 5(r/F„,)2/5 + e. (4) 

The original suggestion of Rosenfeld to use 5 = 3.0 and 
e = 0 [U is shown by the dashed curve. Some im¬ 
provement can be observed when choosing (5 = 3.2 and 
e = —O.I, as documented in Figure [TJ These values are 
therefore adopted throughout this paper. Note that al¬ 
though the functional form Q with the exponent | pro¬ 
vides reasonable accuracy for Yukawa fluids in a wide 
range of k, it can be not the best choice for each single 
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FIG. 1: (Color online) Thermal component of the reduced 
excess energy of Yukawa fluids versus the reduced coupling 
parameter F/Fm. Symbols correspond to the numerical simu¬ 
lations for different values of the screening parameter k [3, [13. 
The (red) solid curve is the fit of Eq. Q with <5 = 3.2 and 
e = —0.1. The (black) dashed line corresponds to the same 
functional form, but with 5 = 3.0 and e = 0, as suggested 


value of K. For instance, in the case of one-component- 
plasma (limiting case /t = 0 of Yukawa systems), the 
exponent i is known to deliver better accuracy [23 - IHI| . 

Equations ©-(jll) with the proper expression for rm(K) 
(see below) provide a simple and accurate tool to esti¬ 
mate the excess energy of Yukawa fluids. The excess free 
energy can then be obtained by the standard integration 

/.,(«,r) = £ ““(pro jp. (5, 

Note that in case of the non-zero value of the parameter 
6, this integral is diverging logarithmically. A simple con¬ 
ventional procedure to avoid this divergence is to start 
integration from F = I in ([S]) and add the correspond¬ 
ing value /ex(K, I) i- The contribution from the weakly 
coupled region fex{K, 1) is in fact of minor importance 
for F ^ 1. The values of fexi^, 1) have been tabulated 
in Refs. @,[13, interpolation is straightforward. This is 
not done here, because the main interest is the excess 
pressure and its derivative at strong coupling, where the 
contribution from the weak coupling region is negligible. 

Using the dimensionless quantities k and F instead of 
particle density and temperature, the excess pressure can 
be obtained from [IJ, [13] 


Pex(K,r) 


r 5/ex _ K dfex 
3 (?F 3 dn ' 


( 6 ) 


After the straightforward manipulation we can rewrite 
this as 


Pex{K,T) = Pst{K,r) + -Uth(K,r) 


5k duth{K, F) 


6 Ok 


( 7 ) 
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where pst(K,r) = ir[Mf(K) — K.dMi{K)/dK\ is the ex¬ 
cess pressure associated with the static component of the 
internal energy and the two last terms in (l7|) come from 
the thermal component of the internal energy. The static 
contribution can be easily evaluated and yields a compact 
expression 0 




6 [/tcosh(K) — sinh(K)]^ 


( 8 ) 


To evaluate the thermal contribution, the functional de¬ 
pendence rin(K) should be specified. Several approx¬ 
imate methods to locate the fluid-solid phase transi¬ 
tion, based on the properties of the interaction poten¬ 
tial alone, have been discussed in the literature 
Some of them are applicable to a wide range of inter¬ 
actions, including the Yukawa case, and can be used 
for this purpose. Being focused on Yukawa interaction 
here, we adopt a simple expression proposed by Vaulina 
et al. [13, im 




172 exp(a«;) 
aK,+ ’ 


(9) 


where the constant a = (47r/3)^/^ ~ 1.612 is the ra- 
tio of the mean interparticle distance A = rip ' to the 
Wigner-Seitz radius a. Equation ([9|) is in rather good 
agreement with the numerical data in the regime k ^ 5 
addressed in this study. Equations O, 0, (0), and ([9]) 
allow to evaluate the excess pressure and hence the com¬ 
pressibility factor. 


Z(K,r) = 1-|-pex(«:,r), (10) 


of Yukawa fluids for any given pair k and E. This consti¬ 
tutes the main result of this paper. In the Appendix 
the explicit expressions for the compressibility factor and 
isothermal compressibility modulus are provided, which 
are convenient for practical calculations. 


IV. COMPARISON WITH PREVIOUS STUDIES 

To demonstrate the accuracy of the present approach, 
the compressibility factor of a single component Yukawa 
fluid has been evaluated in a wide range of coupling 
strength and compared with the results available in the 
literature. This comparison is shown in Table H The 
first three columns contain the location of the system 
state point in terms of k, T, and E/Tni, respectively. 
The fourth column lists the results from the MC sim¬ 
ulations performed by Meijer and Frenkel (ME) [3] and 
tabulated in Ref. [ll|, which serve as a reference data 
here. The fifth column contains the results obtained 
using the modified Rogers-Young (RY) integral equa¬ 
tion [HI- This method requires less numerical efforts 
compared to the original RY implementation, but the 
thermodynamic consistency is achieved only on a dis¬ 
crete grid of points and for chosen fitting functions [HI- 


TABLE I: Compressibility factor (reduced pressure) of a sin¬ 
gle component Yukawa fluid in a wide range of coupling. The 
first two columns specify the location of the system state point 
on the (k, P) plane. The third column lists the values of the re¬ 
duced coupling strength L/Fm (note that the first point may 
correspond to supercooled liquid). The remaining columns 
contain the values of Z obtained using MC simulations by 
Meijer and Frenkel (MF) [3| (^mf), DRY method by Tejero 
et al. 0 (^dry), present approach (.^present), and its static 
component (Zst). For details see the text. 


K 

r 

T/Tm 

Zmf 

^DRY 

■^present 

Zst 

1.800 

396.9 

1.03 

102.492 

102.751 

102.526 

99.856 

1.860 

383.9 

0.95 

89.606 

89.846 

89.567 

86.902 

1.923 

371.4 

0.87 

78.148 

78.387 

78.145 

75.487 

1.984 

360.0 

0.80 

68.640 

68.865 

68.637 

65.987 

2.049 

348.6 

0.73 

59.889 

60.091 

59.895 

57.256 

2.117 

337.5 

0.66 

52.133 

52.307 

52.150 

49.523 

2.182 

327.3 

0.60 

45.711 

45.862 

45.707 

43.095 

2.238 

319.2 

0.56 

41.041 

41.176 

41.002 

38.404 

2.306 

309.7 

0.51 

35.954 

36.072 

35.903 

33.324 

2.348 

304.2 

0.48 

33.204 

33.314 

33.184 

30.618 

2.398 

297.9 

0.45 

30.294 

30.394 

30.249 

27.699 

2.532 

282.1 

0.37 

23.780 

23.855 

23.741 

21.237 

2.631 

271.5 

0.32 

20.016 

20.069 

19.989 

17.524 

2.778 

257.1 

0.26 

15.705 

15.722 

15.682 

13.279 

3.050 

234.2 

0.18 

10.400 

10.343 

10.418 

8.144 


It is therefore referred to as the discretized RY method 
(DRY). The sixth column summarizes the results ob¬ 
tained using the present approach, Eq. ini. The last 
column provides the static contribution to the compress¬ 
ibility factor, Zst(K,T) = 1 -|-pst(K,r). 

The simple approximation discussed in the present pa¬ 
per is in excellent agreement with the precise results from 
MC simulations by Meijer and Frenkel. The deviations 
from their numerical data do not exceed few parts in 
one thousand in all cases considered and are significantly 
smaller than that at strong coupling, near the fluid-solid 
phase transition. Present approach demonstrates better 
accuracy than the DRY integral equation method. The 
static component of the compressibility provides reason¬ 
able estimate of the actual compressibility near freezing 
(^ 3% underestimation), but becomes progressively less 
accurate when coupling decreases. 

Since the numerical data tabulated in Ref. 0 are lim¬ 
ited to a relatively narrow range of the screening parame¬ 
ter K, we have performed additional comparison with the 
MC simulations of Yukawa systems on a hypersphere, 
performed by Caillol and Gilles (CG) [l3]- The high¬ 
est value of the coupling parameter F = 100 investi¬ 
gated in Ref. 0 has been chosen for detailed compar¬ 
ison. The values of the compressibility factor obtained 
in MC simulation and those calculated with the help of 
the expression 0 are listed in Table El for a number of 
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TABLE II: Compressibility factor (reduced pressure) of the 
coirventional Yukawa fluid (with ireutralizing background) for 
a fixed coupling parameter L = 100 and various screening 
parameters k. The first two columns specify the location of 
the system state point on the (k, V /Tm) plane. The third 
and fourth columns contain the values of Z obtained using 
MC simulations by Caillol and Gilles (CG) and using the 
present approach (present), respectively. 


K 

r / r ,. 

ZcG 

■^present 

0.1 

0.58 

- 28.1399 

- 28.1391 

0.2 

0.58 

- 28.0233 

- 28.0291 

0.4 

0.57 

- 27.5706 

- 27.5819 

0.6 

0.54 

- 26.8387 

- 26.8495 

0.8 

0.50 

- 25.8645 

- 25.8698 

1.0 

0.45 

- 24.6883 

- 24.6885 

1.4 

0.35 

- 21.9057 

- 21.9122 

2.0 

0.22 

- 17.3432 

- 17.3674 

2.5 

0.14 

- 13.7956 

- 13.7996 

3.0 

0.08 

- 10.8036 

- 10.7379 

3.5 

0.05 

- 8.43152 

- 8.26628 

4.0 

0.03 

- 6.60858 

- 6.34905 

5.0 

0.01 

- 4.15975 

- 3.80722 


K-values (0.1 < k < 5.0). To keep the original nota¬ 
tion of Ref. we have added the contribution of the 
neutralizing medium to the compressibility factor (see 
Appendix |B|, that is why the compressibility is negative. 
The agreement is again excellent at sufficiently strong 
coupling. The relative deviations do not exceed a tiny 
fraction of a percent as long as k ^ 3.0, where T/Tm has 
dropped well below 0.1. For higher k, deviations increase 
and reach several percent. This apparently indicates that 
the contribution from the weak coupling region has to 
be properly accounted for in this regime. Overall, the 


present approach demonstrates excellent performance at 
least in the regime k ^ 5 and T/Tm ^0.1 and therefore 
can hnd application in many practical situations relevant 
to colloidal systems and complex plasmas. 


V. CONCLUSION 

Simple approach to estimate the internal energy, pres¬ 
sure, and compressibility modulus of three-dimensional 
Yukawa fluids across coupling regimes has been put for¬ 
ward. Explicit analytical expressions for these quantities 
have been derived (see Appendix lAI) . which demonstrate 
excellent agreement with precise results from MC simu¬ 
lations at strong coupling. These expressions are directly 
applicable to single component Yukawa systems, modifi¬ 
cations to describe charged particles in the neutralizing 
medium are trivial (see Appendix |B|) . The obtained re¬ 
sults can be particularly useful in connection with study¬ 
ing wave phenomena in strongly coupled complex (dusty) 
plasmas, since simple and accurate equation of state is 
often required for such studies. For this reason, the main 
focus here has been on pressure-related quantities, al¬ 
though other thermodynamic functions can apparently 
be also estimated using the proposed approach. The ap¬ 
proach is also likely to be relevant for other simple fluids 
with soft repulsive interactions, when the static compo¬ 
nent of the internal energy dominates over the thermal 
one. This and related issues are left for future work. 
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Appendix A: Practical expressions for compressibility factor and isothermal compressibility modulus 


To keep some generality we consider the coefficients 5 and e, entering Eq. (|T]) for the thermal component of the 
excess energy, as unspecified for a moment. The explicit expression for the compressibility factor is 


Z{K,r) — ^1 -I- -b 


Fk^ 

6 [KCOsh(K) — sinh(K)]^ 


I 

3 



2/5 


(Al) 


where 


fz{x) 


+ x'^ +2x + 2 
x'^ + 2x + 2 


(A2) 


The isothermal compressibility modulus is related to the compressibility factor via fi = Z + {T/3)(dZ/dr) — 
(k/3)(9Z/9k). This yields 


/i(K,r) = (^1 -b 


rK®sinh(K) r ^ f 1 

9 [k cosh(K) — sinh(K)]^ 45 vrmy ^ 


(A3) 
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where 




+ 14x® + SSx"' + 76x^ + 136x^ + 136a; + 68 
(a;2 + 2a; + 2)2 ' 


(A4) 


The dependence rni(K) is given by Eq. Regarding the coefficients S and e, we suggest to use 6 = 3.2 and e = — 0.1. 


Appendix B: Effect of the neutralizing medium 


The contributions from particle-particle correlations 
and particle-background interactions are additive. In 
particular, the excess energy associated with the pres¬ 
ence of neutralizing medium is [3^ 


3r kT 

~ ~ T’ 


(Bl) 


The first term represents the excess (free) energy of the 
medium that, on average, neutralizes the charge of the 
particles while the second term gives the (free) energy 


of the sheath around each particle. Note that this latter 
term does not contribute to the excess pressure, as can be 
clearly seen from Eq. (jH]). It, therefore, has no effect on 
the compressibility modulus, too. To account for the ef¬ 
fect of the neutralizing medium, one simply needs to add 
the term —3r/2 k^ to the compressibility factor, Eq. (lAIL 
and the term —3r /k^ to the isothermal compressibility 
modulus, Eq. (jA3l) . Note that the contribution of the 
neutralizing medium is negative and dominant at strong 
coupling, implying that the excess energy, pressure, and 
compressibility are also negative in this regime. 
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